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A COMPARISON BETWEEN THE STELLAR AND DYNAMICAL MASSES OF SIX GLOBULAR CLUSTERS 
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ABSTRACT 

We present the results of a comprehensive analysis of the structure and kinematics of six Galactic globular 
clusters. By comparing the results of the most extensive photometric and kinematical surveys available to date 
with suited dynamical models, we determine the stellar and dynamical masses of these stellar systems taking 
into account for the effect of mass segregation, anisotropy and unresolved binaries. We show that the stellar 
masses of these clusters are on average smaller than those predicted by canonical integrated stellar evolution 
models because of the shallower slope of their mass functions. The derived stellar masses are found to be also 
systematically smaller than the dynamical masses by ^ 40%, although the presence of systematics affecting 
our estimates cannot be excluded. If confirmed, this evidence can be linked to an increased fraction of retained 
dark remnants or to the presence of a modest amount of dark matter. 

Subject headings: globular clusters: general — Methods: data analysis — Stars: kinematics and dynamics — 
Stars: luminosity function, mass function — Stars: Population II 
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L INTRODUCTION 

Among the stellar system zoo there is a widely recognized 
discrepancy between the photon-emitting mass (constituted 
by stars and gas) and dynamical mass (including the contri- 
bution of dark matter). From galaxy clusters to dwarf galax- 
ies the mass of the baryonic component, estimated through 
the conversion of light in mass, has been found to be signif- 
icantly smaller than that estimated via kinematical consider- 
ations (see Tollerud et al. 2011 for a recent review). The 
commonly accepted solution to this problem invokes the pres- 
ence of a large amount of dark, non-baryonic matter driving 
the kinematics of these systems without contributing to their 
luminosity. Within this framework, globular clusters (GCs) 
are considered a notable exception: although they are the 
next step down from the smallest stellar systems containing 
dark matter (the Ultra Compact Dwarf galaxies; Mieske et al. 
2008) and brighter than the most dark matter dominated sys- 
tems (the Ultra Faint Galaxies; Simon et al. 2011) they are 
not believed to contain dark matter (Baumgardt et al. 2009). 
This conclusion is mainly justified by the evidence that dy- 
namical mass-to-light (M/L) ratios of GCs are ^25% smaller 
than those predicted by simple stellar population models that 
assume a canonical Initial Mass Function (IMF; McLaughlin 
& van der Mai-el 2005; Strader et al. 2009, 2011). However, 
these models do not account for dynamical effects such as the 
preferential loss of low-mass stars due to energy equiparti- 
tion which alters the shape of the present-day mass function 
(PDMF) and consequently the stellar M/L ratio. In fact, all 
the past determinations of PDMFs in GCs derived slopes shal- 
lower than that of a canonical (e.g. Kroupa 2001) IMF (Piotto 
& Zoccali 1999; Paust et al. 2010; De Marchi et al. 2010). In 
a recent paper Kruijssen & Mieske (2009) modeled the evo- 
lution of the mass function of 24 GCs including the effects of 
dynamical dissolution, low-mass star depletion, stellar evolu- 
tion and stellar remnants finding a substantial agreement be- 
tween the prediction of their model and the dynamical M/L 
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ratio estimated through the cluster kinematics. They conclude 
that dynamical effects are likely responsible for the observed 
discrepancy between stellar and dynamical masses. The flood 
of available data from deep photometric and spectroscopic 
surveys make now possible to test this issue observationally 
at an higher level, sampling the mass function of GCs down 
to the hydrogen burning limit (Paust et al. 2010) and their 
velocity dispersion profile up to many half-mass radii (Lane 
et al. 2010a). This gives the unprecedented opportunity to 
derive stellar masses from direct star counts instead of con- 
verting light in mass, obtaining an estimate which does not 
require the assumption of a M/L ratio. 

In this paper we compare the results of deep Hubble Space 
Telescope (HST) photometric observations and wide-field ra- 
dial velocities available from recent public surveys with a set 
of dynamical models with the aim of deriving the amount of 
stellar and dynamical mass in a sample of six Galactic GCs. 

2. OBSERVATIONAL MATERIAL 

For the present analysis we make use of three different 
databases: 

• The set of publicy available deep ACS@HST photo- 
metric catalogs of the "globular cluster treasury project" 
(Sarajedini et al. 20073 

• The radial velocities obtained by the AAOmega@AAT 
survey of GCs (Lane et al. 201 ifl 

• The surface brightness profiles derived by Trager et al. 
(1995fl 

The photometric database consists of high-resolution HST 
observations of a sample of 65 Galactic GCs. The database 
has been constructed using deep images secured with the 
Advanced Camera for Surveys (ACS) Wide Field Channel 
through the F606W and F814W filters. The field of view 
of the camera (202" x 202") is centered on the clusters cen- 
ter with a dithering pattern to cover the gap between the 
two chips, allowing a full coverage of the core of all the six 

http://www.astro.ufl.edu/ata/public_hstgc/ 
' http://cdsarc.u-strasbg.fr/viz-bin/qcat 7J/A+A/530/A3 1 
ftp://cdsarc.u-strasbg.fr/pub/cats/J/AJ/109/218 
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GCs considered in our analysis. This survey provides deep 
color-magnitude diagrams (CMD) reaching the faint Main Se- 
quence (MS) of the target clusters down to the hydrogen burn- 
ing limit (atMv < 10.7) with a signal-to-noise ratio S/N > 10. 
The results of artificial star experiments are also available 
to allow an accurate estimate of the completeness level and 
photometric errors. A detailed description of the photomet- 
ric reduction, astrometry and artificial star experiments can 
be found in Anderson et al. (2008). Since all the consid- 
ered GCs have relaxation times shorter than their ages (Mc 
Laughlin & van der Marel 2005), the effects of mass seg- 
regation are expected to be important. It is therefore useful 
to constraint the radial variation of the MF outside the half- 
mass radius. Unfortunately, the ACS treasury survey cov- 
ers only the central part of the cluster at distances < 2' from 
the cluster centers. Additional deep data sampling the exter- 
nal part of the cluster are available only for NGC288. They 
consist of a set of ACS images observed under the program 
GO-12193 which are centered at ^ 5' away from the clus- 
ter center In particular, 3 x 200 H- 1 x 15 i long exposures 
have been taken through the F606W filter and 3 x 150 i H- 
1 X 10 i long exposures through the F814W one. All im- 
ages were passed through the standard ACSAVFC reduction 
pipeline. Data reduction has been performed on the individ- 
ual pre-reduced (.fit) images using the SExtractor photometric 
package (Bertin & Amouts 1996). For each star we measured 
the flux contained within a radius of 0. 1 arcsec (correspond- 
ing to 2 pixel ^FWHM) from the star center We preformed 
the source detection on the stack of all images while the pho- 
tometric analysis has been performed independently on each 
image. Only stars detected in two out of three long-exposures 
or in the short ones have been included in the final catalogue. 
We used the most isolated and brightest stars in the field to 
link the aperture magnitudes at 0.5 arcsec to the instrumen- 
tal ones, after normalizing for exposure time. Instrumental 
magnitudes have been transformed into the VEGAMAG sys- 
tem by using the photometric zero-points by Sirianni et al. 
(2005). Finally, each ACS pointing has been corrected for 
geometric distortion using the prescriptions by Hack & Cox 
(2001). We preformed artificial star experiments on the sci- 
ence frames following the prescriptions described in Sollima 
et al. (2007). In summary, a set of artificial stars have been 
simulated using the Tiny Tim model of the ACS PSF (Krist 
et al. 2010) and added to all images at random positions 
within 36 x 36px cells centered on a grid of 100 x 50 knots 
along the x and y directions of the two ACS chips (a single 
star for each cell). The F814W magnitude of each star has 
been extracted from a luminosity function simulated adopt- 
ing a Kroupa (2001) mass function which has been converted 
in F814W magnitude using the mass-luminosity relation of 
Dotter et al. (2007). The corresponding F606W magnitude 
has been derived by interpolating along the cluster mean ridge 
Une of the F8 14W-(F606W-F8 14W) CMD. We performed the 
photometric reduction on the simulated frames with the same 
procedure adopted for the science frames producing a catalog 
of ~ 10^ artificial stars. 

The radial velocity database has been derived by Lane et 
al. (2009, 2010a, 2010b, 2011) from a large number of 
spectra of Red Giant Branch (RGB) stars observed in fields 
centered on 10 Galactic GCs. Spectra were obtained with 
the multi-fiber AAOmega spectrograph mounted at the Anglo 
Australian Telescope with the 1700D and the 1500V gratings 
on the red and blue arm, respectively. With this configuration 



spectra covering the Ca II triplet region (8340-8840 A) and the 
swathe of iron and magnesium lines around ^ 5200A were 
obtained with a resolution of R= 10000 and R=3700 for the 
red and blue arms, respectively. A detailed description of 
the reduction procedure and radial velocities estimates can 
be found in Lane et al. (2009, 2010a, 2010b, 2011). For 
the present work we adopted the radial velocities extracted 
with the RAVE pipeline since they provide a better estimate 
of the radial velocity uncertainty when compared with the 
available high-resolution spectroscopic studies (see Bellazz- 
ini et al. 2012). In the analysis we used only stars satisfying 
the same membership criterion adopted by these authors and 
within d <2 rit from the cluster center (see Sect. S); they are 
132 in NGC288, 164 in NGC5024, 190 in NGC6218, 357 in 
NGC6752, 676 in NGC6809 and 193 in NGC7099. 

The surface brightness profiles were taken from the Trager 
et al. (1995) database. They were constructed from gener- 
ally inhomogeneous data based mainly on the Berkeley Glob- 
ular Cluster Survey (Djorgovski & King 1984). The sur- 
face brightness profile of each cluster has been derived by 
matching several set of data obtained with different tech- 
niques (aperture photometry on CCD images and photo- 
graphic plates, photoelectric observations, star counts, etc.). 
Although these inhomogeneities may represent a drawback 
of the analysis (see Sect. 16.1b . this database represents the 
unique collection of wide field surface brightness profiles cov- 
ering the entire extent of the target clusters up to their tidal 
radius. 

The target clusters analyzed in this work were selected 
among the 10 objects in common between the three above 
databases requiring that /) the photometric observations reach 
the magnitude of stars with masses M ^ 0. 15 Mq with a com- 
pleteness level > 50%, and //) the systemic rotation does not 
exceed the 10% of the cluster velocity dispersion. This last 
requirement is necessary since rotation affects the determina- 
tion of the actual velocity dispersion (and consequently the 
dynamical mass) both providing a source of rotational kinetic 
energy and introducing an azimuthal variation of star counts 
and kinematics (Wilson 1975) which will not be included in 
our models. Six clusters satisfy the above criteria, namely: 
NGC288, NGC5024, NGC6218, NGC6752, NGC6809 and 
NGC7099. 

3. DYNAMICAL MODELS 

Following the prescriptions of Gunn & Griffin (1979), the 
phase-space distribution of stars is given by the contribution 
of eight groups of masses with the distribution function 
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where E and L are, respectively, the energy and angular mo- 
mentum per unit mass, v,. and v, are the radial and tangential 
components of the velocity, the effective potential i/i is the dif- 
ference between the cluster potential at a given radius r and 
the potential at the cluster tidal radius tp = cj)- (j),, Ai and ki 
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where po is the central cluster density and 
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is the core radius (King 1966). The potential at each radius is 
determined by the Poisson equation 



(3) 



Equations|2]and[3]have been integrated once assumed a value 
of the adimensional potential at the center Wo outward till the 
radius r, at which both density and potential vanish. The total 
mass of the cluster M is then given by the sum of the masses 
of all the groups 



M = ^ niiNi 



n 2 8 



Figure 1. Selection boxes adopted for the population of single stars (from N\ 
loNi) and binaries (At,,,) of NGC 288. The F814W-(F606W-F814W) color- 
magnitude diagram is overplotted. The magnitude corresponding to the 50% 
completeness is marked by dashed line. 



are scale factors for each mass group, and ctk is a normaliza- 
tion term. For each cluster we defined the eight mass bins (all 
covering equal mass intervals at different ranges) to sample 
the entire mass range from the hydrogen burning limit to the 
mass at the end of the Asymptotic Giant Branch evolution set 
by comparison between the CMD and the isochrones by Dot- 
ter et al. (2007)(see Fig. [T]i. The dependence on mass of the 
coefficients A,- determines the degree of mass segregation of 
the cluster. Since the half-mass relaxation times of all the six 
GCs of our sample are significantly smaller than their ages 
(McLaughlin & van der Marel 2005) we adopted A{i) oc m,- 
(Michie 1963). The parameter ra determines the radius at 
which orbits become more radially biased. As oo, mod- 
els become isotropic. As usual, the above distribution func- 
tion has been integrated to obtain the number density and the 
radial and tangential components of the velocity dispersion of 
each mass group: 
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The above equations can be written in terms of dimensionless 
quantities by substituting 
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As a last step, the above profiles have been projected on the 
plane of the sky to obtain the surface mass density 
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and the line-of-sight velocity dispersion 
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In the above models the shape of the density and velocity dis- 
persion profiles are completely determined by the parameters 
(Wo, fa,Ni) while their normalization is set by the pair of pa- 
rameters (rc,M). 

For each choice of the above parameters a large number 
(N,o, = 10*) of particles has been randomly extracted from 
a continuous mass function which reproduces the number 
counts in each mass bin and distributed across the cluster ac- 
cording to their masses by interpolating through the density 
profiles of the various mass groups. For each particle the tan- 
gential and radial components of its velocity have been then 
extracted from the distribution function defined by eq. [Tjand 
projected in an arbitrary direction Q. The masses of the par- 
ticles have been converted in F814W and Johnson V mag- 
nitude adopting the mass-luminosity relation of Dotter et al. 
(2007). Isochrones has been selected adopting the metallic- 
ities by Carretta et al. (2009aj3 and their magnitudes have 
been converted from the absolute to the apparent plane adopt- 
ing the distances and reddening by Paust et al. (2010fl while 
ages have been chosen to match the turn-off morphology of 
each cluster. We adopted the reddening coefficients listed by 
Sirianni et al. (2005). The F606W magnitude has been instead 
assigned by interpolating on the MS mean ridge line of each 

^ this is equivalent to randomize the direction of the line of sight 
^ The metallicity of NGC 5024, not included in the Carretta et al. (2009a) 
sample, has been taken from Arellano Ferro et al. (201 1). 

' The distance and reddening of NGC 5024 and NGC 6218 are not in- 
cluded in the Paust et al. (2010) sample, and have been taken from the latest 
version of the Harris catalog (Harris 1996, 2010 edition). 
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cluster in the F814W-(F606W-F814W) CMD. The same pro- 
cedure has been performed for a fraction ft, of binaries and 
f renin of dark remnants whose masses and magnitudes have 
been assigned following the prescriptions described in Sect. 
13.11 A synthetic Horizontal Branch (HB) has been simulated 
for each cluster using the tracks by Dotter et al. (2007), tuning 
the mean mass and mass dispersion along the HB to repro- 
duce the observed HB morphology. Photometric errors and 
completeness corrections have been included in the follow- 
ing way: to each simulated particle we associated an artificial 
star within 10" and with input F814W magnitude within 0.25 
mag. If the articial star is recovered its corresponding shifts 
in color and F814W magnitude have been then applied to the 
particle, while the particle has been removed if the associated 
artificial star has not been recovered in the output artificial star 
catalog. To account for Galactic field interlopers, we used the 
Galaxy model of Robin et al. (2003). A catalogue covering 
an area of 1 square degree around each cluster center has been 
retrieved. A subsample of stars has been randomly extracted 
from the entire catalogue scaled to the ACS field of view and 
the V and I Johnson-Cousin magnitudes were converted into 
the ACS photometric system by means of the transformations 
of Sirianni et al. (2005). The correction for photometric er- 
rors and incompleteness have been applied to this sample as 
for cluster stars. So, at the end of the above procedure ac- 
cording to a given choice of input parameters a mock cluster 
with masses, projected velocities and magnitudes accounting 
for observational effects has been produced. 

3.1. Binaries and Dark Remnants 

The presence of binaries can influence both the luminos- 
ity and the mass estimate of a cluster Indeed, while the two 
components of the binary system contribute to the total bud- 
get of light and mass in the same way as if they were single 
stars, their mutual interaction has two important effects: /) 
since the two components are not resolved the binary system 
can share the same portion of the CMD with single stars of 
different mass (Malkov & Zinnecker2001), and //) from a dy- 
namical point of view they behave as a single massive object 
so that, as energy equipartition is established, they present on 
average lower velocities and a more concentrated distribution 
with respect to single stars. The above characteristics make 
them more resistant to the process of mass-loss via evapo- 
ration and tidal shocks. It is therefore important to account 
for the presence of a fraction of unresolved binaries to prop- 
erly model the PDMF of the cluster To model the population 
of binaries we extracted Ntin = 2fbNioi stars from the Kroupa 
(2001) IMF in the mass range 0.1 < M/Mq < 7 and paired 
them randomly. From the resulting sample we extracted a 
subsample of pairs according to a "chance of pairing" calcu- 
lated by imposing that /) the distribution of mass ratios f(q) 
(where q = m2/m\ is the ratio between the masses of the sec- 
ondary and the primary components of each binary system) in 
the mass range 1 < M /Mq < 7 results the same observed by 
Fisher et al. (2005) in their sample of spectroscopic binaries in 
the solar neighborhoocQ and //) the distribution of systemic 
masses must be the same of single stars in the common mass 
range. The accepted pairs were added to a library of binaries 
while the rejected ones are "broken" and their components 

Although the evolution of the characteristics of binaries in the solar 
neighborhood is expected to be very different from that in GCs we adopted 
this sample since it provides the only available constraint on the binaries mass 
ratios distribution. 



have been used for the next iteration. The procedure has been 
repeated until all stars are included in the library. Then, all 
binaries whose primary component has a mass exceeding the 
maximum allowed for an evolving cluster star have been re- 
moved. The location and velocity of binaries within the clus- 
ter has been assigned according to their systemic mass (see 
Sect. O. The relative projected velocity of the primary com- 
ponent has been added to the projected velocity of each binary 
according to the formula 



liraisini , ,„ . , 
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McConnachie & Cote (2010) 



where a\ is the semimajor axis of the primary component, P 
is the orbital period, e is the eccentricity, / is the inclination 
angle to the line-of-sight, 6 is the phase from the periastron, 
and uj is the longitude of the periastron. We followed the pre- 
scriptions of Duquennoy & Mayor (1991) for the distribution 
of periods and eccentricities. The semimajor axis has been 
calculated using Kepler's third law: 
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where nii and m2 are the masses of the primary and sec- 
ondary components. We removed all those binaries whose cor- 
responding semimajor axes lie outside the range «,„,„ < ai + 
fl2 < flflMA where fl„„„ is linked to the radius of the secondary 
component (according to Lee & Nelson 1988) and a„,ax is set 
physically by the maximum separation beyond which the bi- 
nary becomes unbound due to stellar interactions within the 
cluster (see eq. 1 of Hills 1984). The distribution of the angles 
(i, 9, uj) has been chosen according to the corresponding prob- 
ability distributions (Prob(0 cx sin/;Prob(6') cx 6'"';Prob(a;) = 
constant). The magnitudes of the two components have been 
then estimated as for single stars and the overall magnitude 
has been calculated from the sum of their fluxes in the dif- 
ferent passbands. The photometric errors and completeness 
correction have been then applied following the same proce- 
dure adopted for single stars (see Sect. |3). 

The number of dark remnants (white dwarfs, neutron stars 
and black holes) has been estimated adopting for the precur- 
sors a Kroupa (2001) IMF normalized to the number of stars 
in the most massive stellar birP^. We then adopted the initial- 
final mass relation and the retention fractions (for the initial 
velocity kicks) by Kruijssen (2009). We applied an addi- 
tional correction to account for the selective removal of low- 
mass remnant as a result of tidal effects by imposing the same 
PDMF in the mass range in common with single stars. 

4. METHOD 

The choice of the best model that reproduces the structure 
and kinematics of each cluster has been made by comparing 
simultaneously four observables measured on the databases 
and in the simulated mock cluster: 

• The surface brightness profile; 

• The F814W luminosity function; 

' ' Such a choice is driven by the fact that massive stars are the most resis- 
tant against tidal losses which are expected to decrease the relative fraction of 
low-mass stars. On the other hand, stellar evolutionary effects are expected to 
be negligible since the majority of the stars in this mass bin are located below 
the turn-off point. 
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Figure 2. Comparison between three observables of NGC 288 (filled points) and the corresponding model prediction (solid lines). Upper-left panel: surface 
brightness profile; upper-right panel: velocity dispersion profile; bottom panels: F814W luminosity function in the inner (left) and outer (right) cluster regions. 



• The fraction of binaries; 

• The velocity dispersion profile. 

The surface brightness profile of the mock cluster has been 
derived by summing the Johnson V fluxes of particles con- 
tained within circular concentric annuli of variable width and 
dividing the resulting flux by the annulus area. The widths 
of the annuli have been chosen to contain at least 20 stars to 
minimize Poisson noise while maintaining a good resolution 
in distance. A of the fit has been then calculated. 

The comparison of the F814W luminosity function and 
of the binary fraction have been performed simultaneously 
by comparing the number counts {N"''^) in nine regions of 
the F814W-(F606W-F814W) CMD defined as foUows: eight 
F814W magnitude intervals have been defined in correspon- 
dence to the mass bins adopted in the cluster modeling (see 
Sect. [3]) and including all stars with colors within 3 times 
the photometric error corresponding to their magnitudes. A 
binary region has been defined within four boundaries: the 
selected stars are those contained in magnitude between the 



loci of binaries with primary star mass wi = 0.45 Mq (faint 
boundary) and mi = 0.75 Mq (bright boundary), and in color 
between the MS ridge line (blue boundary) and the equal- 
mass binary sequence (red boundary) both red-shifted by 3 
times the photometric error (see Fig. [T] see also Sollima 
et al. 2007). We excluded the innermost 0.5' of the core- 
coUapsed clusters NGC6752 and NGC7099 as in these ex- 
tremely crowded regions the completeness level drops below 
50% even at relatively bright magnitudes. 

The velocity dispersion profile of the cluster has been com- 
pared using the likelihood merit function 

where r,, v, and Sv, are the distance to the cluster center, 
the radial velocity and the corresponding uncertainty of the 
i-th observed star, v is the mean velocity of the sample and 
(Ty{ri) is the velocity dispersion predicted for RGB stars by the 
model at r, . We excluded the outermost region of the cluster 
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at li > 2 r/, to minimize the effect of tidal heating that can 
inflate the velocity dispersion in the outskirts of our clusters 
(Johnston et al. 1999; Kupper et al. 2010). 

We adopted an iterative algorithm that starts from an ini- 
tial guess of the parameters (WQ,rc,M,ra,Ni,fi,), generates the 
mock cluster and calculates the observables described above. 
Corrections to the parameters are then calculated using the 
method described below and the procedure is repeated until 
convergence. In each iteration the following steps are per- 
formed sequentially: 



1. 



The value of Wo and rc providing the best fit to the 
surface brightness profile are searched by mini- 
mizing the x~ statistic. At the end of this step a 
new mock cluster is generated leaving the other 
parameters unchanged; 

The relative fraction of stars in the eight mass bins 
A^, is adjusted by multiplying them for corrective 
terms which are proportional to the ratio between 
the number of objects in each bin of the observed 
and of the simulated CMDs. 



N;=Ni 



mock 



obs 



where A^"*' and A^;'""^* are the number of stars 
counted in the i-th bin in the observed CMD and 
in simulated one, N'j is the updated value of A^, and 
is a softening parameter, set to 0.5, used to avoid 
divergence. An updated value for the binary frac- 
tion fh is also set using the same method. Then, a 
new mock cluster is generated and the procedure 
is repeated from point (1) until convergence. 

• The above steps are repeated for different values of fa 
and the merit function / is calculated by comparison 
with the radial velocities to determine the best-fit value 

of fa. 

The procedure usually converges after 10-20 iterations inde- 
pendently on the initial guess of parameters. Convergence is 
set when the various parameters start to fluctuate around the 
best-fit value with an amplitude (typically < 1%) linked to the 
degree of degeneracy among the various solutions. The un- 
certainty in the derived parameters has been calculated as the 
standard deviations of such fluctuations quadratically com- 
bined with the Poisson noise of star counts (only for the case 
of the M,Ni and ft, parameters). The mass of the model can 
be constrained in two independent ways: /) by best fitting the 

total number counts (so imposing Ym=\ N"""^'' = X^Li -^f''') or 
ii) by matching the amplitude of the velocity dispersion. The 
former estimate gives the stellar mass (M*), the latter the dy- 
namical one (Mdvn)- 

In Fig. |2]the comparison between the considered observ- 
ables and the prediction of the best-fit model are shown for 
the case of NGC288. It is clear that the model well-reproduces 
the surface brightness and velocity dispersion profiles as well 
as the F814W luminosity functions measured in the two ra- 
dial samples considered for this cluster In this cluster, where 
ancillary observations are available in an external region, the 
model that best fit the inner sample observables also fits both 
the shape and the normalization of the luminosity function in 
the external sample. Although this agreement cannot guaran- 
tee the full adequacy of models and observables in the whole 




log M,,yMj, 

Figure 3. Derived fraction of binaries as a function of tlie dynamical mass. 

sample of clusters, it represents a sanity check supporting the 
validity of the assumptions made on mass segregation and on 
the star counts-surface brightness conversion. 

5. RESULTS 

The procedure described above allows to determine a set of 
parameters for the target clusters, in particular: their structural 
parameters, degree of anisotropy, binary fraction, PDMF and 
their stellar and dynamical masses. The best-fit parameters 
for each cluster are listed in Table 1 . For completeness, the 
M/Ly ratios, calculated from the V magnitudes by van den 
Bergh et al. (1991), are also listed. The structural parameter 
derived here are in good agreement with those estimated by 
previous studies (e.g. McLaughlin & van der Marel 2005). 
We note that all the clusters of our sample are best fitted by 
isotropic (or marginally anisotropic) models. However, the 
relatively large uncertainties in the binned velocity dispersion 
measures do not allow a meaningful analysis of this issue: the 
likelihood merit function / (see Sect. IHl vary by less than 10% 
over a relatively large range of fa values in all cases. 

The main results on the binary fraction, PDMF and mass 
are discussed in the following sections. 

5.1. Binary fractions 

The binary fractions listed in Table 1 are remarkably lower 
than those estimated in previous studies for the clusters in 
common (Sollima et al. 2007; Milone et al. 2012). This 
is related to the distribution of mass ratios adopted in this 
work that is different from those used in previous studies. The 
above difference can be well responsible for the observed dis- 
crepancies since in binaries with low mass ratios the contri- 
bution of the secondary component to the total flux is negli- 
gible, and they appear in the CMD almost indistinguishable 
from single stars. Therefore, the adoption of a mass ratio 
distribution skewed toward low (high) mass ratios produce 
a larger (smaller) correction for hidden binaries. To verify 
the consistency of our analysis with the previous works we 
calculated the fraction of binaries in NGC288 adopting a flat 
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Table 1 

Best-fit parameters of the six target clusters 



Name 


Wo 


rc 






h 


fcore 
Jb 


Q 




log M, 




M,/Lv 








pc 




% 


% 




Mq 


Mq 


Mq/Lv.q 


Mq/Lv/,0 


NGC 288 


6.2 


1.90 


5.40 


11 


3.3±0.2 


4.8±0.3 


-1.04±0.05 


4.83±0.07 


4.73±0.I 


1.3±0.4 


1.0±0.4 


NGC 5024 


11 


0.55 


2.99 


+00 


4.3±3.0 


11.5±7.2 


-1.08±0.07 


5.60±0.07 


5.61±0.1 


1.5±0.3 


1.5±0.4 


NGC 6218 


7.7 


1.05 


1.50 


12 


3.2±0.2 


5.8±0.4 


-0.22±0.02 


4.90±0.05 


4.74±0.1 


1.0±0.3 


0.7±0.2 


NGC 6752 


13 


0.42 


0.51 


+00 


0.6±0.1 


3.6±0.6 


-1.41±0.09 


5.45±0.04 


5.20±0.1 


2.4±0.5 


1.3±0.4 


NGC 6809 


5.5 


2.50 


3.32 


4 


2.5±0.2 


3.4±0.3 


-0.67±0.01 


5.19±0.03 


5.08±0.1 


1.6±0.2 


1.3±0.2 


NGC 7099 


13.5 


0.13 


0.34 


+00 


1.7±0.1 


12.0±0.9 


-1.14±0.04 


5.28±0.05 


5.08±0.1 


2.1±0.4 


1.3±0.3 



T r 




NGC6218 

\ , , , \ , 1 , \ , , , \ , , , I 1 , 

-0.8 -0.6 -0.4 -0.2 

log M/Mg 



Figure 4. The present-day mass function of the six target clusters are shown. 
A vertical shift has been added for each cluster for clarity. 

mass-ratios distribution over the entire mass range instead of 
the assumptions described in Sect. 13.11 In this case the over- 
all binary fraction turns out tobe6.5±0.3% corresponding to 
9.6 ± 0.4% in the cluster core, in agreement within the errors 
with the estimates by Sollima et al. (2007) and Milone et al. 
(2012). It is worth noting that both the PDMF and the mass 
estimates are almost unaffected by this change: in this case in 
fact the slope of the PDMF varies of less than A a ~ -0.07 
and both the dynamical and luminous masses remain the same 
within AlogM/MQ < 0.01. So, although the adopted mass- 
ratio distribution is somehow arbitrary, it influences only the 
absolute fraction of binaries having a negligible effect on the 
relative cluster-to-cluster ranking, on the luminous and dy- 
namical masses and on the shape of the PDMF. 

In Fig. [3] the fraction of binaries of our sample of clus- 
ters are plotted as a function of the derived cluster dynamical 
mass. Once excluded NGC 5024 which has a large uncer- 
tainty in the derived binary fraction, a clear anticorrelation is 
noticeable with massive clusters hosting a smaller fraction of 
binaries. This result has been already observed in Milone et 
al. (2008, 2012) and SolUma et al. (2010) on a larger sample 
of clusters. 

5.2. Present-Day Mass Function 

In Fig. |4]the global PDMFs of the six clusters are shown. 
The PDMFs can be well represented by single power laws 
whose slopes are all comprised within the range -0.2 < a < 
-1.45 being significantly less steep than the IMF of Salpeter 



(1955; a = -2.35) and, in aU but the case of NGC6752, of 
Kroupa (2001; a = -1.3 at M < 0.5 Mq with a mean slope 
of < a >= -1.7 between 0.3 < M/Mq < 0.8). The exteme 
value for NGC 6218 (a = -0.22 ± 0.02) confirms the pecu- 
liar flatness of the PDMF of this cluster already noticed by 
De Marchi et al. (2006). For the two clusters NGC 288 
and NGC7099 in common with Paust et al. (2010) the de- 
rived slopes agree within the combined uncertainties of the 
two works (^ 0.2). In the same way, the slopes derived for 
NGC 6809 and NGC 7099 agree to those estimated for these 
clusters by Piotto & Zoccah (1999) and Piotto et al. (1997), 
respectively. We find no signs of a turnover at low masses in 
all cases, although a convex shape of the PDMF is visible in 
NGC 5024. It is however worth noting that the mass at which 
the PDMF is expected to culminate is ^ 0.3 Mq (Paresce & 
De Marchi 2000), so this last conclusion depends only on the 
relative fraction of stars in the lowest-mass bin. However, the 
estimate in this bin is by far the most uncertain since i) the 
mass-luminosity relation in this mass range is largely uncer- 
tain, ii) the completeness in this bin is low and uncertain, and 
Hi) some contamination from faint unresolved galaxies falling 
in the large selection box of this mass bin (see Fig. [U could be 
present. So, a firm conclusion on the presence of a turnover 
in the PDMF in these clusters cannot be established. 

6. DYNAMICAL VERSUS STELLAR MASS 

The main result of this analysis is the measure of the dy- 
namical and stellar masses of the target clusters. Previous 
estimates of the dynamical mass of GCs have been made by 
Mandushev et al. (1991), Pryor & Meylan (1993), Lane et al. 
(2010a) and Zocchi et al. (2012) which have respectively 3, 5, 
6 and 3 clusters in common with the present study. The mean 
differences and standard deviations between the dynamical 
masses estimated by these authors and those estimated here 
are A log{M/MQ) =0.13(cr =0.15), -0.03(0.17), 0.10(0.26) 
and -0.09(0.14^3 respectively, indicating a good agreement 
within the errors without significant systematics. In the right 
panel of Fig. |5]the derived stellar masses are compared with 
those reported by McLaughlin & van der Marel (2005) for the 
four clusters in common using the M/Ly values predicted by 
the population synthesis model of Bruzual & Chariot (2003) 
adopting a Kroupa (2001) IMF. The stellar masses derived 
here are in reasonable agreement within the errors with those 
derived by these authors when considering individual clusters, 
although on average our estimates are smaller by ^ 20%. This 
is expected since the Kroupa (2001) IMF is steeper than the 
PDMFs shown in Fig. |4] overestimating the number of low- 
mass stars that contribute significantly to the clusters' mass 
budget and only marginally to their light. 

For the Zocchi et al. (2012) estimates we adopted those deiived 
with the models, considered more reliable by these authors. The 
mean(dispersion) of the difference in log (M/Mq) with respect to the es- 
timates made adopting the King (1966) models is 0.16(0.23). 
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Figure 5. Left panel: comparison between the stellar and dynamical masses of the six target clusters. Right panel: comparison between the stellar masses 
measured in the present work and in McLaughlin & van der Marel (2005). The one-to-one relation is marked by the dashed line in both panels. 



From an inspection of Table 1 it is suddenly evident a dis- 
crepancy between the dynamical and stellar masses of five out 
of six clusters (with the exception of NGC 5024) with the dy- 
namical masses systematically larger than the stellar ones by 
^40% . This is shown in the left panel of Fig. |5] where the 
two estimates are compared. 

The observed difference between the dynamical and stellar 
masses can be also visualized in Fig. |6]and|7]where the veloc- 
ity dispersion profile and of the F814W luminosity function 
of the six target clusters are compared with the prediction of 
the best-fit models with the appropriate M* and M^,,,. It is 
clear that, apart from the case of NGC5024, models are not 
able to fit simultaneously both the velocity dispersion profile 
and of the F814W luminosity function. Such a discrepancy 
is significant from a statistical point of view: a test indi- 
cates that the probability that the above difference occurs by 
chance in the entire sample of clusters is smaller than 3%. 
Such a discrepancy can be due to either observational biases 
in the determination of masses or to a physical reason. In the 
next sections we will analyze both possibilities. 

6.1. Possible observational biases 

To understand if the observed deficiency of stellar mass is 
real or due to an effect of the adopted technique we must 
search for some systematics in our analysis leading to an over- 
estimate of the dynamical masses or to an underestimate of the 
stellar ones in all the six clusters of our sample. 

The dynamical mass mainly follows from the fit of the ve- 
locity dispersion profile, therefore any systematics affecting 
the velocity dispersion can alter its determination. An ob- 
vious possibility is that the velocity errors reported by Lane 
et al. (2011) are underestimated. In fact, according to eq. 
m the observed velocity spread is given by the convolution 
of the "true" velocity dispersion and the observational errors. 
So, for a given observed velocity spread, the larger are the 
errors the smaller is the estimated "true" velocity dispersion 
and, consequently, the best-fit dynamical mass is expected to 



be smaller than estimated. However, the comparison with a 
high-resolution sample shown by Bellazzini et al. (2012) in- 
dicates a scatter compatible with (and even smaller than) the 
combined uncertainties of the two datasets. This is also in 
agreement with what is reported by Lane et al. (2011) who 
claim that the velocity uncertainties of their sample could be 
overestimated by ^ 20%, making the difference even more 
striking. 

The stellar mass depends mainly on the normalization be- 
tween the observed and predicted F814W luminosity func- 
tion and on the correction for the fraction of mass not cov- 
ered by the ACS field of view. In all the analyzed clus- 
ters the samples of stars involved in the fit consist of more 
than 15000 stars, ~90% of them ai-e above the 80% level of 
completeness. So, the Poisson error due to number counts 
and to the error on the completeness correction are negligi- 
ble (AlogM^/MQ < 0.01). A good correction for the radial 
coverage is supposed to be provided by the excellent fit of the 
surface brightness profile. However, we must stress that such 
profiles are constituted by a compilation of inhomogeneous 
data, often merging aperture photometric measurements (trac- 
ing the distribution of the brightest and most massive RGB 
stars) in the central part with star counts of (low-mass) MS 
stars in the outskirts. As stars with different masses have also 
different radial profiles (as a result of mass segregation), this 
introduces a large source of uncertainty in the fit. In this con- 
text, the simultaneous agreement between the predicted and 
observed F8 14W luminosity function of NGC 288 in both the 
internal and external fields is encouraging. 

6.2. Effect of assumptions 

If the discrepancy is not linked to any observational bias, 
it can be related to a physical effect not (properly) consid- 
ered in our models. There are many possible mechanisms 
that may affect the mass determination. In this Section we 
evaluate the impact of the various input taking the cluster 
NGC288 as tester. In this cluster the observed mass discrep- 
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Figure 6. Comparison of the observed velocity dispersion profile (top panels) and the F814W luminosity function (bottom panels) of the clusters NGC288 (left), 
NGC5024 (middle) and NGC6218 (right) with the coiTesponding models with the best-fit M^yn (solid lines; red in the onHne version) and M, (dashed lines; blue 
in the online version). The velocity dispersion measures not included in the analysis (at d > 2 r/,) are marked with open points in upper panels. 



ancy (AlogM/MQ = 0.1) is smaller than the average of the 
six analysed GCs (AlogM/MQ = 0.18) so that care should be 
taken in the general interpretation of the derived results. 

As discussed in the previous section, the estimate of the dy- 
namical mass is sensible to all the parameters affecting the 
velocity dispersion. Radial anisotropy can in principle distort 
the velocity dispersion profile, resulting in a different normal- 
ization of mass. As discussed in Sect. |5] we take into ac- 
count for this effect in our models, but nevertheless the rela- 
tively large uncertainties in the velocity estimates do not al- 
low a firm constraint on this quantity. To test this effect, 
we repeated the fit imposing the maximum degree of radial 
anisotropy that still ensures stability (see Nipoti et al. 2002). 
The largest effect on the derived dynamical mass has been 
found to be A log Mdyn/MQ < 0.03. Binaries can also inflate 
the velocity dispersion because of the presence of the addi- 
tional orbital motion of the binary components. Also in this 
case, while our code accounts for such an effect, a different 
assumption of the distribution of mass ratios might produce a 



difference. However, the test performed in Sect. 15.1 [ indicates 
that while the above assumption can affect significantly the 
fraction of binaries, the resulting mass is only marginally af- 
fected (A log Mdy„/MQ < 0.01). Another important effect is 
the heating due to the tidal interaction of the clusters with the 
Milky Way: the periodic shocks occurring every disk crossing 
and pericentric passage transfer kinetic energy to the cluster 
stars inflating its velocity dispersion, particularly in the out- 
skirts (Gnedin & Ostriker 1997). While this effect is limited 
in our case since we excluded by the fit the outermost por- 
tion of the velocity dispersion profile, Kupper et al. (2010) 
showed that significant heating can occur also within the half- 
mass radius of lose clusters. To test this hypothesis, we run 
an N-body simulation to follow the structural and kinematic al 
evolution of a cluster with the orbital and structural character- 
istics of NGC 288. This is the less dense cluster of our sam- 
ple and should therefore be among those where tidal shocks 
are expected to produce the largest effect. We used the last 
version of the coUisionless N-body code gyrfalcON (Dehnen 
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Figure 7. Same of Fig. |6]for NGC6752 (left), NGC6809 (middle) and NGC7099 (right). 



2000). The cluster has been represented by 61,440 particles 
distributed according to the best-fit multimass model of this 
cluster (see Table 1). The mean masses of th e p articles have 
been scaled to match the observed cluster masS- We adopted 
a leapfrog scheme with a time step of Af = 1.45 x 10'*yr and 
a softening length of 0.43 pc (following the prescription of 
Dehnen & Read 2011). Such a relatively large time step and 
softening length do not affect the accuracy of the simulation 
because we are interested in the tidal effects which occur on 
a timescale significantly shorter than the relaxation time (the 
timescale at which collisions become important). The clus- 
ter was launched within the three-component (bulge -i- disk 
-I- halo) static Galactic potential of Johnston et al. (1995) 
with the energy and angular momentum listed by Dinescu 
et al. (1999) and its evolution has been followed for 0.852 
Gyr (corresponding to 4 orbital periods). The projected den- 
sity and velocity dispersion profiles at the end of the simu- 

The adopted scaling does not affect the results of the simulation as the 
tidal effects we are interested to study depends on the overall cluster potential 
and not on the individual star masses. 



lation are shown in Fig. [8] It is clear that the final veloc- 
ity dispersion deviates from the prediction of the steady-state 
model at a distance of ^ 4' ~ r/, to the center By apply- 
ing the same technique described in Sect. |4]to the outcome 
of the simulation we overestimate the mass of the cluster by 
A/ogM^„,/M0~O.O6. 

The stellar mass can be underestimated because of an im- 
proper conversion between F814W magnitude and mass or 
because of the presence of non-luminous matter Considering 
that the mass discrepancy is of the order of ^40%, uncertain- 
ties on distance, reddening and age can be excluded as possi- 
ble causes: to reproduce the above mass difference one should 
adopt for these clusters ages <4 Gyr, distance moduli shifted 
by ^2.5 mag or reddening differences of A£'(B-y) ^0.2, 
many orders of magnitudes larger than their statistical uncer- 
tainties. Uncertainties in the mass-luminosity relation can be 
also responsible for part of the discrepancy. In this regard, 
it is worth noting that while different stellar evolution mod- 
els all agree for masses M > 0.1 Mq (Alexander et al. 1997; 
Baraffe et al. 1997; Dotter et al. 2007, Marigo et al. 2008) 
significant differences can be present in the least massive bin. 
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Figure 8. Projected density (left panel) and velocity dispersion (right panels) profiles calculated from the last snapshot of the N-body simulation. The quantities 
at the end of the simulation are marked with filled dots. The best-fit static model is oveiplotted in both panels with dashed lines. 



Also in this case, however, the effect of such uncertainty is 
not expected to exceed few percent in the estimated mass. A 
larger effect is instead produc ed by the assumption on dark 
remnants. As described in Sect. 13.11 we followed the prescrip- 
tion by Kruijssen (2009) for the retention fraction of these 
objects and a Kroupa (2001) IMF for their precursors. To 
estimate the extent of the impact of this assumption we con- 
sidered the extreme case of a 100% retention fraction for all 
the remnant. In this ideal case, the stellar mass would increase 
by A log M^,/Mq ^ 0.08. The adoption of a different slope 
for the IMF has an even larger effect: the entire mass dis- 
crepancy could be completely removed by assuming a slope 
of a = -1 .5 (instead of -2.3) for masses M > O.8M0, because 
the shallower the slope of the IMF is, the larger is the relative 
fraction of massive stars which terminated their evolution and 
are now in form of dark remnants. 

7. DISCUSSION 

By comparing the most extensive photometric and spectro- 
scopic surveys with suitable dynamical models we derived the 
global binary fraction, the PDMF of six Galactic GCs as well 
as their stellar and dynamical masses. The approach adopted 
here has the advantage to not require an a priori assumption 
of the M/L ratio to estimate the stellar mass (as usually done 
in previous works; e.g. McLaughlin & van der Marel 2005) 
which is instead estimated via direct star counts in the CMD. 

We confirm the anticorrelation between binary fraction and 
cluster mass already found in previous studies (Milone et al. 
2008, 2012; Sollima et al. 2010) and explained as the con- 
sequence of the increased efficiency of the binary ionization 
process in massive clusters. Indeed, in massive clusters both 
the number of collisions and the mean kinetic energy of stellar 
encounters is larger resulting in a higher probability of ioniza- 
tion of multiple systems (Sollima 2008). 

The PDMFs of the six clusters are well-represented by sin- 
gle power laws, although we cannot exclude a turnover at 
masses M < 0.2 Mq, because of the uncertainty in the rel- 



ative fraction of the least massive stars. 

The stellar masses derived for these clusters are on average 
smaller than those predicted by population synthesis models. 
This is a consequence of the fact that the IMF adopted by 
these models is steeper that the actual PDMF of the consid- 
ered clusters. This has been previously predicted by Kruijssen 
& Mieske (2009) as the result of the preferential loss of low- 
mass stars during the cluster evolution. 

We found a discrepancy between the stellar and dynamical 
masses of five clusters. In particular, dynamical masses are on 
average ^40% larger than stellar ones. Such a discrepancy is 
statistically significant when considering the entire sample of 
six clusters and could not be detected by previous studies be- 
cause of the overestimate of the stellar mass mentioned above. 
Unfortunately, we cannot exclude that such a discrepancy is 
due to the presence of systematics affecting our estimates. In 
particular, a possible source of systematics is constituted by 
the inhomogeneity of the surface brightness profiles sample 
of Trager et al. (1995). However, if the detected difference 
would be real, a number of physical interpretations could be 
given. After a careful test on the effect of various assump- 
tions and physical processes we found that a significant spu- 
rious increase of the estimated dynamical masses is given by 
the tidal heating which can reduce (but not eliminate) the ob- 
served discrepancy. On the other hand, the assumption on the 
distribution and retention of dark remnant has the largest im- 
pact on the estimated stellar masses. Part (up to 75%) of the 
observed mass difference can be indeed explained assuming 
that most the remnants, including the most massive neutron 
stars and black holes, are retained. This could be the case if 
the cluster mass at the moment of the SNe 11 explosion was 
large enough to trap the remnants, whose kinetic energy has 
been increased by the velocity kick following the explosion. 
The mass necessary to retain the majority of massive rem- 
nants is of ^ 10^ Mq (see Fig. 1 of Kruijssen 2009). Such 
a large mass is comparable with the initial mass predicted by 
D'Ercole et al. (2008) to ignite the self-enrichment process 
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observed in most GCs (see Carretta et al. 2009b). More- 
over, massive remants are expected to be more resistant to 
tidal stripping since tend to sink toward in the cluster core 
as a result of mass segregation. It is however worth noting 
that a retention fraction larger than 80% would require initial 
masses and an early mass loss history quite extreme, which 
has indeed been proposed in theoretical works to be the result 
of severe tidal shocking in the primordial environment of GCs 
(Elmegreen 2010; Kruijssen et al. 2012). An even larger ef- 
fect is provided by the slope of the high-mass end of the IMF 
of the precursors which can drastically increase the amount 
of non-luminous mass. In this regard, a recent analysis by 
Marks et al. (2012) suggests that the slope of the IMF in the 
range M > M3 could be actually shallower than the canonical 
Salpeter (1955) and Kroupa (2001) values with a dependence 
on density and metallicity. In this case, a slope a ~ -1 .5 (in- 
stead of q; = -2.3 by Kroupa 2001) could account for the en- 
tire difference between the estimated stellar and dynamical 
masses. Such a variation would have strong implications on 
the evolution of the GC MF: in fact, while the largest mass 
contribution would be due to the increased number of white 
dwarfs, the consistent increase of the fraction of massive rem- 
nants (black holes and neutron stars) makes these objects the 
main drivers of the MF evolution. It is not clear if such IMF 
could actually evolve in the PDMF observed in these clusters. 

Another possibility is that GCs contain a modest fraction 
of non-baryonic dark matter. In this case, it is possible that 
what we see now is the remnant of a larger halo lost during 
the interaction with the Milky Way (see Saitoh et al. 2006). 
Some amount of dark matter is also expected if some GCs are 
the remnants of past accretion events, as previously suggested 
by Freeman & Bland-Hawthorn (2002). However, given the 
involved uncertainties, all these hypotheses are merely spec- 
ulative. Further studies will be needed to verify the above 
possibilities. 
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